BARRETO et al. 2022: Habitat loss and biodiversity gain (BIOCON-S-22-00746)
Model formulas
For each of our response variables we run three candidate models for each scale:
Model 0 <- ~1, absence of effect
Model 1 <- percentage of forest cover at 3 and 5km radius
Model1.2 <- non-linear forest cover at 3 and 5km radius (quadratic model)
Loading data and gathering into one datasets calculated in script ‘data_wrang_GIT.Rmd’, content: 1) “beetles”: raw data for species occurrence per landscape/site; 2) “hab.data”: data containing classes of habitat association; 3) “env.data”: environmental variables, forest cover percent at 3 and 5km radius buffer 4) “div_hab”: diversity (alpha, beta and gamma) between classification habitat association 5) “rc.df”: Beta RC calculations 6) “ab_hab”: mean and median abundance between groups of habitat association 7) “st.hab_ab”: abundance data between classification habitat association per site
| Landscape | gamma_FS | gamma_NFS | alpha_FS | alpha_NFS | betad_FS | betad_NFS | ab.land_FS | ab.land_NFS | brc.fsNA | brc.nfsNA | brc.abund.fsNA | brc.abund.nfsNA | fc_3km | fc_5km |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 148 | 11 | 9 | 3.75 | 2.25 | 7.25 | 6.75 | 66 | 40 | -0.2978571 | 0.1920000 | -0.0130130 | 0.4388388 | 38.52523 | 43.86337 |
| 215 | 19 | 7 | 5.62 | 2.00 | 13.38 | 5.00 | 174 | 51 | 0.0928571 | -0.3038095 | 0.2775990 | 0.3995901 | 30.15485 | 25.23160 |
| 263 | 16 | 13 | 6.00 | 4.50 | 10.00 | 8.50 | 393 | 214 | -0.5305000 | -0.0057143 | 0.6195243 | 0.3910577 | 10.10641 | 13.41180 |
| 266 | 11 | 8 | 3.25 | 1.38 | 7.75 | 6.62 | 70 | 12 | -0.2342857 | 0.1323810 | 0.1122074 | 0.1577292 | 30.87525 | 31.16090 |
| 282 | 18 | 10 | 6.62 | 4.00 | 11.38 | 6.00 | 220 | 211 | -0.2835714 | -0.1857143 | 0.2198985 | 0.6894990 | 15.26270 | 16.35308 |
| 291 | 9 | 3 | 3.71 | 1.71 | 5.29 | 1.29 | 150 | 49 | -0.3942857 | -0.5880952 | 0.0863244 | -0.4210878 | 45.74294 | 45.13213 |
| 317 | 12 | 4 | 5.00 | 1.38 | 7.00 | 2.62 | 152 | 56 | -0.3942857 | -0.5371429 | 0.1745674 | -0.4296439 | 48.75805 | 47.85138 |
| 329 | 21 | 13 | 8.75 | 5.50 | 12.25 | 7.50 | 373 | 258 | -0.0653571 | 0.1196429 | 0.9374374 | 0.4348992 | 14.73791 | 13.62588 |
| 333 | 15 | 14 | 4.62 | 4.50 | 10.38 | 9.50 | 159 | 218 | -0.1078571 | 0.1910714 | 0.6737809 | 0.6187974 | 13.71403 | 14.50960 |
| 335 | 15 | 8 | 5.75 | 2.38 | 9.25 | 5.62 | 140 | 113 | -0.2667857 | -0.3271429 | 0.3398041 | -0.9444206 | 23.87913 | 19.21670 |
| 359 | 17 | 16 | 6.25 | 5.25 | 10.75 | 10.75 | 152 | 130 | -0.1714286 | 0.2753571 | 0.7798870 | 0.8147076 | 18.10208 | 17.62431 |
| 399 | 21 | 13 | 8.00 | 5.00 | 13.00 | 8.00 | 433 | 168 | -0.0861905 | 0.2435714 | 0.9606273 | 0.6988774 | 29.73616 | 24.97851 |
RESULTS
We found 28 species to be forest specialists, those that come from the Atlantic Forest domain and exclusively collected in forested habitats; while 23 were non-forest specialists, species from a broader domain and/or less restrictive to forest habitats (hereafter FS and NFS, respectively).
FOREST SPECIALISTS
Gamma diversity
Modelling
g0 <- glm(gamma_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.fs.aic <- AICctab(g0, g5k, g5kq, g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g5k -28.2 65.4 5.1 0.0 3 0.751
## g3k -30.0 68.9 3.3 3.5 3 0.127
## g5kq -28.1 69.9 5.2 4.6 4 0.076
## g0 -33.3 71.9 0.0 6.6 2 0.028
## g3kq -29.6 72.9 3.7 7.6 4 0.017
Summary top model
summary(g5k)##
## Call:
## glm(formula = gamma_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.24788 -0.14644 -0.02457 0.14244 0.30719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.46708 2.02439 10.604 9.26e-07 ***
## fc_5km -0.23174 0.05885 -3.938 0.00278 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03547038)
##
## Null deviance: 0.80342 on 11 degrees of freedom
## Residual deviance: 0.34475 on 10 degrees of freedom
## AIC: 62.361
##
## Number of Fisher Scoring iterations: 3
g.fs <- g5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(g5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1272, p-value = 0.672
## alternative hypothesis: two.sided
hnp::hnp(g5k)## Gamma model
boot::glm.diag.plots(g5k) ## for glm poisson or neg binomPrediction:
Alpha diversity
Modelling
a0 <- glm(alpha_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.fs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a5k -19.5 48.0 3.0 0.0 3 0.538
## a3k -20.5 50.0 2.0 2.0 3 0.200
## a0 -22.4 50.2 0.0 2.2 2 0.175
## a5kq -19.2 52.1 3.2 4.1 4 0.068
## a3kq -20.5 54.7 2.0 6.7 4 0.019
Summary top model
summary(a5k)##
## Call:
## glm(formula = alpha_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.43968 -0.11810 -0.03523 0.09321 0.36169
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.56549 0.96276 7.858 1.38e-05 ***
## fc_5km -0.07518 0.02853 -2.635 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.06116683)
##
## Null deviance: 1.00974 on 11 degrees of freedom
## Residual deviance: 0.62036 on 10 degrees of freedom
## AIC: 44.989
##
## Number of Fisher Scoring iterations: 4
a.fs <- a5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(a5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1023, p-value = 0.672
## alternative hypothesis: two.sided
hnp::hnp(a5k)## Gamma model
boot::glm.diag.plots(a5k)Prediction:
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_FS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.fs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b5k -22.5 54.0 5.6 0.0 3 0.710
## b5kq -21.8 57.3 6.3 3.3 4 0.135
## b3k -24.4 57.7 3.7 3.8 3 0.108
## b3kq -23.3 60.2 4.8 6.2 4 0.031
## b0 -28.1 61.5 0.0 7.6 2 0.016
Summary top model
summary(b5k)##
## Call:
## glm(formula = betad_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.24292 -0.15104 -0.01746 0.05238 0.31092
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.93904 1.28413 10.855 7.46e-07 ***
## fc_5km -0.15787 0.03686 -4.283 0.0016 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03504058)
##
## Null deviance: 0.83647 on 11 degrees of freedom
## Residual deviance: 0.33078 on 10 degrees of freedom
## AIC: 50.972
##
## Number of Fisher Scoring iterations: 4
b.fs <- b5k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(b5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1273, p-value = 0.624
## alternative hypothesis: two.sided
hnp::hnp(b5k) ## for glm Gamma or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b5k) ## for glm Gamma or neg binomPrediction
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc0 4.7 -4.0 0.0 0.0 2 0.451
## rc3kq 8.1 -2.5 3.4 1.5 4 0.213
## rc5k 5.3 -1.6 0.6 2.4 3 0.136
## rc5kq 7.5 -1.2 2.8 2.8 4 0.112
## rc3k 4.9 -0.7 0.2 3.3 3 0.088
Summary of 2nd best model, as the top is of absence of effect (~1)
summary(rc3kq)##
## Call:
## lm(formula = brc.fsNA ~ fc_3km + I(fc_3km^2), data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1824 -0.1016 -0.0105 0.0858 0.2064
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7047653 0.2339873 -3.012 0.0147 *
## fc_3km 0.0431907 0.0182460 2.367 0.0421 *
## I(fc_3km^2) -0.0007821 0.0003078 -2.541 0.0317 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1423 on 9 degrees of freedom
## Multiple R-squared: 0.4366, Adjusted R-squared: 0.3114
## F-statistic: 3.487 on 2 and 9 DF, p-value: 0.07563
rc1.fs <- rc0# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc3kq)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.82544, p-value = 0.816
## alternative hypothesis: two.sided
hnp::hnp(rc3kq)## Gaussian model (lm object)
Prediction:
Total abundance
Modelling
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_FS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_FS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_FS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_FS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.fs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -70.6 150.1 2.1 0.0 3 0.404
## gab0 -72.7 150.7 0.0 0.5 2 0.311
## gab3k -71.2 151.4 1.5 1.2 3 0.218
## gab5kq -70.4 154.6 2.2 4.5 4 0.044
## gab3kq -71.1 155.8 1.6 5.7 4 0.023
Summary 2nd best model, as top model was the one of absence of effect (~1):
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_FS ~ fc_5km, data = data_div, init.theta = 4.550032026,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6556 -0.9899 -0.2919 0.6493 1.8325
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.93867 0.31886 18.62 <2e-16 ***
## fc_5km -0.02505 0.01108 -2.26 0.0238 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.55) family taken to be 1)
##
## Null deviance: 17.357 on 11 degrees of freedom
## Residual deviance: 12.410 on 10 degrees of freedom
## AIC: 147.13
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.55
## Std. Err.: 1.84
##
## 2 x log-likelihood: -141.134
gab.fs <- gab5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1526, p-value = 0.496
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance FS",
"FC at 5km",
round(gab.fs.aic$AICc[1],2),
round(gab.fs.aic$dAICc[1],1),
gab.fs.aic$df[1],
round(gab.fs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "~ 1 (NULL)",
round(gab.fs.aic$AICc[2],2),
round(gab.fs.aic$dAICc[2],1),
gab.fs.aic$df[2],
round(gab.fs.aic$weight[2],2),
"-"),
c(" ", "FC at 3km",
round(gab.fs.aic$AICc[3],2),
round(gab.fs.aic$dAICc[3],1),
gab.fs.aic$df[3],
round(gab.fs.aic$weight[3],2),
round(gab3k.r2$R2_Nagelkerke,2)))Prediction:
bRaup-Crick (abundance-based)
Modelling:
### LM gamma for regular model selection
rc0 <- lm(brc.abund.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc5k 0.1 8.8 3.9 0.0 3 0.602
## rc3k -0.8 10.7 2.9 1.9 3 0.234
## rc0 -3.8 12.8 0.0 4.1 2 0.079
## rc5kq 0.2 13.3 3.9 4.5 4 0.062
## rc3kq -0.8 15.4 2.9 6.6 4 0.022
Summary top model:
summary(rc5k)##
## Call:
## lm(formula = brc.abund.fsNA ~ fc_5km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38964 -0.18081 -0.01911 0.15650 0.50966
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.910180 0.176628 5.153 0.00043 ***
## fc_5km -0.018384 0.006117 -3.005 0.01322 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2626 on 10 degrees of freedom
## Multiple R-squared: 0.4746, Adjusted R-squared: 0.4221
## F-statistic: 9.033 on 1 and 10 DF, p-value: 0.01322
rc2.fs <- rc5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc5k)## Gaussian model (lm object)
Prediction
NON-FOREST SPECIALISTS
Gamma diversity
Modelling:
### LM alpha for regular model selection
g0 <- glm(gamma_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.nfs.aic <- AICctab(g0,
g5k, g5kq,
g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g3k -26.3 61.6 7.6 0.0 3 0.7974
## g5k -28.4 65.9 5.5 4.2 3 0.0958
## g3kq -26.1 65.9 7.8 4.3 4 0.0950
## g5kq -28.4 70.5 5.5 8.9 4 0.0093
## g0 -33.9 73.2 0.0 11.6 2 0.0025
Summary top model:
summary(g3k)##
## Call:
## glm(formula = gamma_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.38295 -0.23805 -0.05072 0.13490 0.38893
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.39301 1.95020 8.919 4.49e-06 ***
## fc_3km -0.28151 0.04841 -5.815 0.000169 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.07298721)
##
## Null deviance: 2.46666 on 11 degrees of freedom
## Residual deviance: 0.71058 on 10 degrees of freedom
## AIC: 58.646
##
## Number of Fisher Scoring iterations: 5
g.nfs <- g3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(g3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.89592, p-value = 0.936
## alternative hypothesis: two.sided
hnp::hnp(g3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(g3k) ## for glm poisson or neg binomPrediction:
Alpha diversity
Modelling:
### LM alpha for regular model selection
a0 <- glm(alpha_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.nfs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a3k -16.1 41.2 5.7 0.0 3 0.497
## a5k -16.4 41.8 5.4 0.6 3 0.365
## a5kq -15.7 45.1 6.1 3.9 4 0.070
## a3kq -15.9 45.5 5.9 4.3 4 0.057
## a0 -21.8 49.0 0.0 7.8 2 0.010
Summary top model:
summary(a3k)##
## Call:
## glm(formula = alpha_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.66790 -0.14831 0.00666 0.10705 0.54796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.65928 0.79667 7.104 3.28e-05 ***
## fc_3km -0.08855 0.02018 -4.388 0.00136 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1090929)
##
## Null deviance: 2.9678 on 11 degrees of freedom
## Residual deviance: 1.1728 on 10 degrees of freedom
## AIC: 38.181
##
## Number of Fisher Scoring iterations: 5
a.nfs <- a3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(a3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.94044, p-value = 0.88
## alternative hypothesis: two.sided
hnp::hnp(a3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(a3k) ## for glm poisson or neg binomPrediction:
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_NFS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.nfs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b3k -24.4 57.7 5.5 0.0 3 0.734
## b3kq -23.8 61.4 6.0 3.7 4 0.118
## b5k -26.2 61.4 3.7 3.7 3 0.115
## b0 -29.9 65.1 0.0 7.4 2 0.019
## b5kq -26.0 65.6 3.9 7.9 4 0.014
Summary top model:
summary(b3k)##
## Call:
## glm(formula = betad_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.71087 -0.18354 -0.04943 0.17316 0.48985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.78046 1.59008 7.409 2.29e-05 ***
## fc_3km -0.19439 0.03898 -4.987 0.000548 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1092755)
##
## Null deviance: 2.9582 on 11 degrees of freedom
## Residual deviance: 1.2094 on 10 degrees of freedom
## AIC: 54.73
##
## Number of Fisher Scoring iterations: 6
b.nfs <- b3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(b3k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.65169, p-value = 0.624
## alternative hypothesis: two.sided
hnp::hnp(b3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b3k) ## for glm poisson or neg binomPrediction:
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -0.4 9.8 2.1 0.0 3 0.365
## rc0 -2.4 10.2 0.0 0.5 2 0.289
## rc5k -0.9 10.9 1.5 1.1 3 0.206
## rc3kq 0.8 12.1 3.2 2.4 4 0.112
## rc5kq -0.6 14.9 1.8 5.2 4 0.028
Summary of 2nd best model, 1st wqs of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29676 -0.21076 -0.06141 0.23644 0.41254
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.279665 0.187805 1.489 0.1673
## fc_3km -0.012984 0.006398 -2.029 0.0699 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2735 on 10 degrees of freedom
## Multiple R-squared: 0.2917, Adjusted R-squared: 0.2209
## F-statistic: 4.118 on 1 and 10 DF, p-value: 0.06989
rc1.nfs <- rc3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Prediction:
Total abundance
Modelling:
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_NFS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_NFS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_NFS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_NFS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.nfs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -63.9 136.7 5.3 0.0 3 0.402
## gab3k -64.2 137.4 4.9 0.7 3 0.285
## gab5kq -62.0 137.7 7.1 0.9 4 0.253
## gab3kq -63.6 141.0 5.5 4.3 4 0.047
## gab0 -69.1 143.6 0.0 6.8 2 0.013
Summary top model:
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_NFS ~ fc_5km, data = data_div, init.theta = 4.180831522,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8524 -0.4383 0.2265 0.4858 0.8634
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.92551 0.33598 17.64 < 2e-16 ***
## fc_5km -0.04812 0.01180 -4.08 4.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.1808) family taken to be 1)
##
## Null deviance: 29.171 on 11 degrees of freedom
## Residual deviance: 12.765 on 10 degrees of freedom
## AIC: 133.73
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.18
## Std. Err.: 1.78
##
## 2 x log-likelihood: -127.727
gab.nfs <- gab5k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.43345, p-value = 0.192
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab5kq.r2 <- r2(gab5kq)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance NFS",
"FC at 5km",
round(gab.nfs.aic$AICc[1],2),
round(gab.nfs.aic$dAICc[1],1),
gab.nfs.aic$df[1],
round(gab.nfs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "FC at 3km",
round(gab.nfs.aic$AICc[2],2),
round(gab.nfs.aic$dAICc[2],1),
gab.nfs.aic$df[2],
round(gab.nfs.aic$weight[2],2),
round(gab3k.r2$R2_Nagelkerke,2)),
c(" ", "Non-linear FC at 5km",
round(gab.nfs.aic$AICc[3],2),
round(gab.nfs.aic$dAICc[3],1),
gab.nfs.aic$df[3],
round(gab.nfs.aic$weight[3],2),
round(gab5kq.r2$R2_Nagelkerke,2)
))Prediction:
bRaup-Crick (abundance-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.abund.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -7.5 23.9 1.8 0.0 3 0.363
## rc0 -9.3 23.9 0.0 0.0 2 0.362
## rc5k -8.0 25.0 1.3 1.1 3 0.214
## rc3kq -7.3 28.4 2.0 4.4 4 0.039
## rc5kq -7.9 29.5 1.4 5.6 4 0.022
Summary 2nd best model, 1st is of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.abund.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24197 -0.18970 0.05611 0.27708 0.52925
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.81906 0.33888 2.417 0.0363 *
## fc_3km -0.02184 0.01155 -1.892 0.0878 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4935 on 10 degrees of freedom
## Multiple R-squared: 0.2635, Adjusted R-squared: 0.1899
## F-statistic: 3.578 on 1 and 10 DF, p-value: 0.08782
rc2.nfs <- rc3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc3k)
plot(sim)## qu = 0.75, log(sigma) = -3.635142 : outer Newton did not converge fully.
testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Prediction:
AICc-model selection table
| Diversity | Model | AICc | dAICc | df | weight | r2 |
|---|---|---|---|---|---|---|
| Gamma FS | FC at 5km | 65.36 | 0 | 3 | 0.75 | 0.58 |
| Alpha FS | FC at 5km | 47.99 | 0 | 3 | 0.54 | 0.4 |
| FC at 3km | 49.97 | 2 | 3 | 0.2 | 0.29 | |
| Beta FS | FC at 5km | 53.97 | 0 | 3 | 0.71 | 0.61 |
| RCbeta (P/A-based) FS | ~ 1 (NULL) | -3.97 | 0 | 2 | 0.45 |
|
| Non-linear FC at 3km | -2.48 | 1.5 | 4 | 0.21 | 0.31 | |
| Total abundance FS | FC at 5km | 150.13 | 0 | 3 | 0.4 | 0.44 |
| ~ 1 (NULL) | 150.66 | 0.5 | 2 | 0.31 |
|
|
| FC at 3km | 151.37 | 1.2 | 3 | 0.22 | 0.33 | |
| RCbeta (abund-based) FS | FC at 5km | 8.78 | 0 | 3 | 0.6 | 0.42 |
| FC at 3km | 10.67 | 1.9 | 3 | 0.23 | 0.32 | |
| Gamma NFS | FC at 3km | 61.65 | 0 | 3 | 0.8 | 0.73 |
| Alpha NFS | FC at 3km | 41.18 | 0 | 3 | 0.5 | 0.63 |
| FC at 5km | 41.8 | 0.6 | 3 | 0.37 | 0.61 | |
| Beta NFS | FC at 3km | 57.73 | 0 | 3 | 0.73 | 0.62 |
| RCbeta (P/A-based) NFS | FC at 3km | 9.75 | 0 | 3 | 0.37 | 0.22 |
| ~ 1 (NULL) | 10.22 | 0.5 | 2 | 0.29 |
|
|
| FC at 5km | 10.9 | 1.1 | 3 | 0.21 | 0.14 | |
| Total abundance NFS | FC at 5km | 136.73 | 0 | 3 | 0.4 | 0.82 |
| FC at 3km | 137.42 | 0.7 | 3 | 0.28 | 0.79 | |
| Non-linear FC at 5km | 137.66 | 0.9 | 4 | 0.25 | 0.93 | |
| RCbeta (abund-based) NFS | FC at 3km | 23.92 | 0 | 3 | 0.36 | 0.19 |
| ~ 1 (NULL) | 23.92 | 0 | 2 | 0.36 |
|
|
| FC at 5km | 24.98 | 1.06 | 3 | 0.21 | 0.12 |
FIGURES
Gamma hab
### Gamma ####
# r2(g.fs) # 0.55
## prediction g2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(g.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(g.nfs) # 0.73
## prediction g1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds2 <- predict(g.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds2$fit
new2$ic.up <- preds2$fit + preds2$se.fit*1.96
new2$ic.low <- preds2$fit - preds2$se.fit*1.96
## Plot gamma richness ~ fc %
(g.hab <- ggplot(data= new, aes(x= fc_5km, y = pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data= data_div, aes(x=fc_5km, y=gamma_FS), col= "#0072B2") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_point(data= data_div, aes(x=fc_3km, y=gamma_NFS), col= "#D55E00") +
xlab(" ") + ylab("Gamma") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=20.5, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.55"))) +
annotate("text", y=19, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.73"))) + labs(tag="a"))Alpha hab
# r2(a.fs) # 0.40
## prediction a2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(a.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(a.nfs) # 0.63
## prediction a1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(a.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Alpha richness ~ fc_3km
(a.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=alpha_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=alpha_NFS), col= "#D55E00") +
xlab("") + ylab("Alpha") +
scale_color_manual(breaks=c("FS", "NFS"),
values=c('FS'='#0072B2', 'NFS'='#D55E00')) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=14,face="bold"),
legend.position = 'bottom') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=13, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.40")))+
annotate("text", y=11, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.63"))) +
labs(tag="b"))Beta hab
### Beta add ####
# r2(b.fs) # 0.58
## prediction b1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(b.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(b.nfs) # 0.65
## prediction b1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(b.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Beta Diversity ~ fc %
(b.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=betad_FS), col= "#0072B2") +
geom_line(data= new2, aes(x= fc_3km, y= pred),
col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=betad_NFS), col= "#D55E00") +
xlab(" ") + ylab("Beta") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=18, x=40,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.58")))+
annotate("text", y=16, x=38, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.65"))) +
labs(tag="c"))g.hab + a.hab + b.habggsave("figures/plot.hab.jpeg",units="cm", width=20, height=8, dpi=300)Total abundance FS & NFS
## FS
# r2(gab.fs) # 0.44
## prediction gab1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(gab.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(gab.nfs) # 0.82
## prediction gab2
new2 <- new
preds <- predict(gab.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot gamma richness ~ forest cover 5km
(gab.hab <- ggplot(new, aes(x=fc_5km, y=pred)) +
geom_line(linetype= "dashed", col= "#0072B2", aes(alpha= 0.1)) +
# geom_ribbon(aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=ab.land_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_5km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.2, fill= "#D55E00") +
geom_point(data= data_div, aes(x=fc_5km, y=ab.land_NFS), col= "#D55E00") +
scale_x_continuous(limits= c(10,50)) +
xlab(" ") + ylab("Total abundance") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
# scale_x_reverse() +
annotate("text", y=300, x=38, alpha= .5, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.44"))) +
annotate("text", y=270, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.82"))) +
labs(tag="d"))
### Beta RC FS & NFS
## FS
# ## NULL ## r2(rc1.fs)
## prediction rc2.2 and rc2
new <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=30))
preds <- predict(rc1.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(rc1.nfs) # 0.25
## prediction rc2
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot RCbeta abundance ~ forest cover at 3
(rc.hab <- ggplot(data= new, aes(x=fc_3km, y=pred)) +
geom_line(data= new,aes(x=fc_3km, y=pred),
col= "#0072B2", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D33E00", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D33E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.nfsNA), col= "#D33E00") +
xlab("Percent forest cover") + ylab(c(expression(bold("β"["RC"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),
legend.position = 'none') + ylim(-1,1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=43, alpha= .5, colour= "#D33E00",
label=expression(paste(~ R^2 ~ "= 0.25"))) +
# annotate("text", y=0.7, x=43, alpha= .5, colour= "#0072B2",
# label=expression(paste(~ R^2 ~ "= 0.30"))) +
# annotate("text", y=0.9, x=43, colour= "#D33E00",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.37"))) +
# annotate("text", y=0.8, x=43, colour= "#0072B2",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.33"))) +
labs(tag="e"))
Beta RC abund FS & NFS
# r2(rc2.fs) # 0.41
## prediction rc2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(rc2.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(rc2.nfs) # 0.21
## prediction rc1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc2.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## RCbeta abund ~ FC%
(rc.ab.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new,aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=brc.abund.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), linetype= "dashed", col= "#D55E00", alpha=.2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.abund.nfsNA), col= "#D55E00") +
xlab(" ") + ylab(c(expression(bold("β"["RC-abundance"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(-1,1.1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=45, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.41"))) +
annotate("text", y=0.75, x=43, alpha= .5, colour= "#D55E00",
label=expression(paste( ~ R^2 ~ "= 0.21"))) +
labs(tag="f"))Final plot
(plot.hab <- (g.hab + a.hab + b.hab) /
(gab.hab + rc.hab + rc.ab.hab))ggsave("figures/figure2.jpeg",units="cm", width=28, height=15, dpi=300)